Method of determining fluid flow

ABSTRACT

A method of determining fluid flow in a volume containing two or more fluid components, including determining a pressure field for the volume. One or more streamlines are determined from the pressure field, and the fluid composition is solved for the fluid composition along the, or each, streamline. The pressure may also be solved along the, or each, streamline. The step of solving along the, or each, streamline may be performed using a finite difference technique.

The present invention relates to a method of determining fluid flow, andto processing fluid flow data relating to a hydrocarbon reservoir.Typically, when a hydrocarbon reservoir is exploited one or moreboreholes are drilled into the reservoir and liquid hydrocarbons such asoil and/or gas hydrocarbons is/are extracted from the reservoir throughthe boreholes (which are generally known as “producers”). A fluid,normally water, is introduced into the reservoir at one or more pointsdistant from the producers, to displace liquid or gaseous hydrocarbonsand cause them to be expelled from the reservoir through the producers.A point at which water is introduced into the reservoir is generallyknown as an “injector”.

The extraction of liquid and/or gaseous hydrocarbons at the producer(s)and the injection of water at the injector set up a fluid flow patternin the reservoir. The fluid flow pattern is also influenced by gravity.It is desirable to model the fluid flow within the reservoir, in orderto predict, for example, how the yield of liquid and/or gaseoushydrocarbons extracted by a particular producer will vary over theplanned lifetime of the exploitation of the reservoir (which can be manyyears). Furthermore a reservoir will generally contain a number ofdifferent hydrocarbon compounds, and it is desirable to model the yieldof a particular hydrocarbon compound (generally referred to as a“component”) at a producer as a function of time. Thus, it is preferablenot just to calculate the overall hydrocarbon pressure at a producer,but also to calculate the hydrocarbon composition at a producer.

The mathematical equations governing fluid flow are known, so that amathematical model of fluid flow in a given reservoir may beconstructed. However such a mathematical model is generally not capableof being solved exactly, and it is necessary to use numerical techniquesto determine an approximate solution. There are two principal classes ofnumerical techniques—streamline methods and finite differencetechniques. Finite difference techniques have generally been used forsimulation of fluid flow in a reservoir, although streamline methodshave been used for numerical approximation of fluid flow since the 1800sand have been used in petroleum engineering since the 1950s.

Conventional finite difference techniques suffer from two principaldrawbacks:—numerical smearing and loss of computational efficiency formodels with a large number of grid cells. (In order to represent amulti-producer, geologically heterogeneous reservoir, a model having 10⁵to 10⁶ cells is required.) Finite difference methods based on an IMPES(Implicit Pressure Explicit Saturation) approach are limited to using ashort time step length owing to the limiting CFL condition, which limitsthe maximum distance that fluid can flow during a time step to less thanthe length of a grid cell. The need to use a short time step lengthincreases the number of steps required in a simulation. As a result, thetotal processing time required to process a large model can make itimpractical to apply the method.

An alternative approach is to use a fully implicit finite differencetechnique. This can, in principal, use a longer time step than anIMPES-based finite difference technique. However, it requires theinversion of a much larger matrix and this is a particular disadvantagewith a simulation that incorporates the composition of fluid at aproducer, since the existence of a large number of hydrocarboncomponents will make the matrix very large. Furthermore, the equationsgoverning fluid flow are non-linear and this can place an upper limit onthe allowable time step length. As a result, the processing timerequired for a fully implicit finite difference simulation may also makeit impractical to apply the method to a large model, in particular to acompositional model. This problem has hitherto been addressed by usingan adaptive implicit finite difference technique.

Conventional streamline techniques are based on an IMPES method. Inthese techniques the fluid pressure is solved implicitly, andstreamlines are computed based on this solution for the fluid pressure.A streamline is the path that a particle would take through athree-dimensional space (in this case the reservoir) for a pressuresolution at a given time. In this method, the three-dimensional domainof the reservoir is decomposed into many one-dimensional streamlines,and fluid flow calculations are performed each streamline. A streamlinetechnique may have a high time step length, which reduces the processingtime required. However, a conventional streamline technique assumes,despite the movement of fluid that occurs in the reservoir, that thepressure is constant over a time step of fluid movement and thisassumption leads to a lack of connection between the pressure field andthe movement of fluids. This can cause instabilities and limit themaximum time step length that can be used.

In an example of conventional streamline technique, the componentconservation equations can be written in terms of moles as:

$\begin{matrix}\begin{matrix}{{{\frac{\partial}{\partial t}( {\phi\;{\sum\limits_{{j = g},{w\;\rho}}{x_{i,j}b_{j}S_{j}}}} )} + {\overset{\rho}{\nabla}{\cdot {\sum\limits_{{j = o},w,g}{x_{i,j}b_{i}{\overset{\rho}{v}}_{j}}}}}} = q_{i}} & \; & {{i = 1},\ldots\mspace{11mu},n_{c}}\end{matrix} & (1)\end{matrix}$where φ is the porosity of the reservoir material, x_(i,j) is the molefraction of component i in the j^(th) phase, b_(j) is the molar densityfor phase j, S_(j) is the saturation for phase j, and q_(i) is the wellproduction rate of component i. The saturation of a phase is thefraction of the total fluid volume for a given phase so, for example, ifhalf the fluid at a given point is oil then the oil saturation is 0.5.The subscript i, where 1≦i≦n_(c) denotes the particular hydrocarboncomponent, there are n_(c) hydrocarbon components in total in thereservoir fluid, and j denotes the phase (gas, liquid etc). The velocityv_(j) of the j^(th) phase is defined as:

$\begin{matrix}{{\overset{\rho}{v}}_{j} = {- {\lambda_{j}( {{\overset{\mu}{\nabla}P_{j}} - {\rho_{j}\overset{\rho}{g}}} )}}} & (2)\end{matrix}$where λ_(j) is the mobility for phase j, P_(j) is the pressure of thephase j, and ρ_(j) is the mass density of phase j.

The streamline concept is based on an operator splitting technique forthe conservation equations, which separates an equation for the pressurefrom the conservation equations. Assuming that capillary forces arenegligible so that P_(c)=0, then a pressure equation may be written as

$\begin{matrix}{{{{- \overset{\rho}{\nabla}} \cdot {\overset{\rho}{v}}_{t}} + Q} = {{\phi\; C_{t}\frac{\partial p}{\partial t}} + {( {\sum\limits_{{j = o},w,g}{{\overset{\rho}{v}}_{j}c_{j}}} ) \cdot {\nabla p}}}} & (3)\end{matrix}$where Q is the total well production rate which equals the total wellinjection rate, v_(t) is the total velocity of fluid, C_(t) is the fluidtotal compressibility; c_(j) is the compressibility of phase j and p isthe pressure.

Again assuming that the capillary pressure P_(c)=0 (in which case theP_(o)=P_(w)=P_(g)), then the total velocity is given by

$\begin{matrix}{{\overset{\rho}{v}}_{t} = {{{- \lambda_{t}}{\overset{\mu}{\nabla}P_{o}}} + {( {\sum\limits_{j}{\lambda_{j}\rho_{j}}} )\overset{\rho}{g}}}} & (4)\end{matrix}$where P_(o) is the oil phase pressure, λ_(t) is the total mobility,λ_(j) is the mobility of phase j, ρ_(j) is the mass density of phase jand g is the acceleration due to gravity.

Equations (3) and (4) are used to find the pressure distribution for agiven component distribution. The pressure equation is solved implicitlyusing a finite difference method to give a pressure field. From thispressure field, cell fluxes are found using equation (4). These fluxesdefine the streamlines. These conventional streamline calculations arewell-known and are described in the literature.

Neglecting the gravitational term in the above equations, the IMPES-typeassumption that the pressure is constant throughout a time step meansthat pressure-derivative terms can be ignored when solving for thesaturation so that, for the concentration equation along a streamline,it is assumed that:

$\begin{matrix}\begin{matrix}{{\overset{\rho}{\nabla}{\cdot {\overset{\rho}{v}}_{t}}} = 0} & \; & {\frac{\partial\phi}{\partial t} = 0}\end{matrix} & (5)\end{matrix}$

The total number of moles per unit volume of component i, m_(i), isdefined as

$\begin{matrix}{m_{i} = {\sum\limits_{{j = o},w,g}{x_{i,j}b_{j}S_{j}}}} & ( {6a} )\end{matrix}$and the convective flow of component i F_(i) is defined as

$\begin{matrix}{F_{i} = {\sum\limits_{{j = o},w,g}{f_{j}x_{i,j}b_{j}}}} & ( {6b} )\end{matrix}$where f_(j) is the fractional flow function of component j and b_(j) isthe molar density of phase j.

If the effect of gravity is neglected, then equation (1) may bere-written as a 1-D equation along the streamlines:

$\begin{matrix}\begin{matrix}{{\frac{\partial m_{i}}{\partial t} + {\frac{1}{\phi}{{\overset{\rho}{v}}_{t} \cdot {\overset{\rho}{\nabla}F_{i}}}}} = q_{i}} & \; & {{i = 1},\ldots\mspace{11mu},n_{c}}\end{matrix} & (7)\end{matrix}$where q_(i) is the injection/production rate of the component i.

The velocity v_(t) defines the streamlines. On the assumption that thepressure in the reservoir remains constant over a time step, then thestreamlines will not change their path during a time step.

In order to account for the effect of gravity, we define the term forthe flow due to gravity by

$\begin{matrix}\begin{matrix}{{{\overset{\mu}{G}}_{i} = {\sum\limits_{{j = o},w,g}{x_{i,j}b_{j}\Gamma_{j}\overset{\rho}{g}}}},} & {\Gamma_{j} = {{\lambda_{j}\rho_{j}} - {f_{j}{\sum\limits_{{k = o},w,g}{\lambda_{k}\rho_{k}}}}}}\end{matrix} & (8)\end{matrix}$which arises from writing v_(t)=v_(t)+Γ_(j)

for each phase j. In equation (8), G_(i) is the flow due to gravitysegregation for component i, Γ_(j) is the gravity segregation strengthfor phase j and the index k denotes oil, water or gas (and is used sincethe usual index j is used earlier in the equation).

Equation (1) may now be re-written to separate the gravity flow termfrom the convective flow term, to give

$\begin{matrix}\begin{matrix}{{\frac{\partial m_{i}}{\partial t} + {\frac{1}{\phi}( {{{\overset{\rho}{v}}_{t} \cdot {\overset{\rho}{\nabla}F_{i}}} + {\overset{\rho}{\nabla}{\cdot {\overset{\rho}{G}}_{i}}}} )}} = q_{i}} & \; & {{i = 1},\ldots\mspace{11mu},n_{c}}\end{matrix} & (9)\end{matrix}$

In equation (9), the ∇.F_(i) term is a 1-D term along the streamlinesdefined by v_(t) and the ∇.G_(i) term is a 1-D term along the gravitylines defined by the gravity vector g. It is possible to split the twoterms, to give:

$\begin{matrix}{{{\frac{\partial m_{i}}{\partial t} + {\frac{1}{\phi}{{\overset{\rho}{v}}_{i} \cdot {\nabla F_{i}}}}} = q_{i}}\mspace{11mu}} & (10) \\{{\frac{\partial m_{i}}{\partial t} + {\frac{1}{\phi}{\overset{\rho}{\nabla}{\cdot {\overset{\rho}{G}}_{i}}}}} = q_{1}} & (11)\end{matrix}$

Equation (10) models the convective flow in the absence of gravity, andmay be solved as described above. Equation (11) models the gravitysegregation of the fluids. The gravity vector in equation (8) isconstant and known, so that the gravity lines are always defined. Thefull 3-D movement of fluids in the reservoir is modelled by solvingequations (10) and (11) sequentially. FIG. 1 illustrates schematicallyhow the saturation is modelled initially along a streamline, and thenalong a gravity line.

The principal steps of a conventional streamline method are shown inFIG. 2. At step 1 a fluid flow model of the reservoir is initialised,and at step 2 an initial estimate of the pressure field in the reservoiris made by solving implicitly for the pressure. At step 3 the initialstreamlines are calculated from the initial estimate of the pressurefield. The saturation is then mapped onto these initial streamlines atstep 4, and the saturation is solved over a time step ΔT along thestreamlines at step 5, using equation (10). Several techniques for thisare known.

The saturation is them mapped from the streamlines onto the gravitylines, at step 6, and the saturation is solved over the time step ΔTalong the gravity lines at step 7, using equation (11). The saturationis then mapped onto grid cells at step 8, so that the saturation at eachgrid point within the reservoir can be output.

To repeat the calculations over a subsequent time step, the time isupdated at step 9, and the pressure field within the reservoir isre-computed from the newly-calculated values for the saturation. Step 3of estimating the streamlines from the pressure is then repeated. Asuitable counter is updated at step 10, and the newly-calculated valuesfor the saturation are then mapped onto the new streamlines. Steps 5 to8 of computing the saturation along the (new) streamlines, mapping thesaturation onto the gravity lines, computing the saturation along thegravity lines and mapping the saturation onto the grid cells are thenrepeated. Steps 9, 10 and 2–8 may then be repeated until the fluid flowhas been simulated over a desired time period.

One advantage of a streamline technique is that the problem ofsimulating the fluid flow in the reservoir has been reduced to a set ofseparate, smaller 1-D problems. As noted above, however, conventionalstreamline techniques assume that the pressure field is constant over atime step (which means that the streamlines are fixed over a time step).This assumption limits the length of the time step that can be used,since effects that depend on pressure changes are ignored. This problemis particularly serious if the pressure within the reservoir is changingrapidly since, in this case, large errors may occur unless the timesteps are be kept short.

A first aspect of the present invention provides a method of determiningfluid flow in a volume containing two or more fluid components, themethod comprising: determining a pressure field for the volume;determining one or more streamlines from the pressure field; and solvingfor the fluid composition along the or each streamline.

A second aspect of the invention provides a method of processing datarelating to a hydrocarbon reservoir containing two or more fluidcomponents, the method comprising the steps of: determining a pressurefield for a volume indicative of the reservoir; determining one or morestreamlines from the pressure field; and solving for the fluidcomposition of the reservoir along the or each streamline.

The present invention thus makes it possible to perform a fullcompositional calculation along a streamline.

In a preferred embodiment the method further comprises solving for thepressure along the or each streamline.

The step of solving the fluid composition along the or each streamlineand, optionally, the step of solving the pressure along the or eachstreamline may be repeated by updating the pressure field, re-computingthe streamlines for the updated pressure field, and solving along there-computed streamlines. This enables the fluid composition, andoptionally the pressure, to be modelled as a function of time.

Where the invention is applied to a hydrocarbon reservoir, it providesinformation about the expected yield at a producer of the reservoir.This may be used to control the production of the reservoir, for exampleto determine where to locate a producer and/or an injector.

In preferred embodiment, the step of solving along the or eachstreamline comprises using a finite difference technique. Thiseliminates the requirement of a conventional streamline technique thatthe pressure is assumed to the constant during a time step. In turn,this allows a longer time step length, and so reduces the number ofsteps required to simulate over a given time period.

A third aspect of the present invention provides an apparatus fordetermining fluid flow in a volume containing two or more fluidcomponents, the apparatus comprising: means for determining a pressurefield for the volume; means for determining one or more streamlines fromthe pressure field; and means for solving for the fluid compositionalong the or each streamline.

A fourth aspect of the present invention provides an apparatus forprocessing data relating to a hydrocarbon reservoir containing two ormore fluid components, the apparatus comprising: means for determining apressure field for a volume indicative of the reservoir; means fordetermining one or more streamlines from the pressure field; and meansfor solving for the fluid composition of the reservoir along the or eachstreamline.

In a preferred embodiment, the apparatus comprised a programmable dataprocessor.

A fifth aspect of the present invention provides a storage mediumcontaining a program for a data processor of an apparatus as definedabove.

Preferred embodiments of the present invention will now be described byway of illustrative example with reference to the accompanying Figures,in which:

FIG. 1 is a schematic illustration of a conventional streamlinetechnique;

FIG. 2 is a flow chart for a conventional streamline technique;

FIG. 3( a) is a flow chart illustrating the principal steps of a methodof the invention;

FIG. 3( b) is a flow chart illustrating one embodiment of the invention;

FIG. 4 is a schematic illustration of one step of a method according tothe present invention;

FIG. 5 illustrates gas production rates predicted by a method of theinvention and a conventional technique for a first reservoir model;

FIG. 6 illustrates oil, gas and water production rates predicted by amethod of the invention and a conventional technique for a secondreservoir model;

FIG. 7 shows a saturation plot for a method of the invention for thesecond reservoir model;

FIG. 8 shows a streamline plot for a method of the invention for thesecond reservoir model;

FIG. 9 illustrates oil production rates predicted by a method of theinvention and a conventional technique for a third reservoir model

FIG. 10 shows a saturation plot for a method of the invention for athird reservoir model;

FIG. 11 shows a streamline plot for a method of the invention for thethird reservoir model;

FIG. 12 shows a streamline plot for a method of the invention for afourth reservoir model;

FIG. 13 is a comparison of the CPU time and RAM capacity required by amethod of the invention and a conventional technique;

FIG. 14 is a block schematic diagram of an apparatus constituting afurther embodiment of the invention; and

FIG. 15 is a block schematic diagram of a programmable dataprocessor-based apparatus constituting a further embodiment of theinvention.

In a streamline calculation of the present invention, the pressure fieldfor a reservoir is calculated at for each time step, and the streamlinesare estimated from this pressure field as in the prior art. However, theone-dimensional solutions along the streamlines are found a fullyimplicit technique or an adaptive implicit technique. The pressure andthe fluid composition are solved for together along each streamline.Physical effects that depend on the changes in pressure can be takenaccount of in the solution along the streamlines, and it is notnecessary to assume that the pressure remains constant over a time step.The present invention is thus able to use long time steps, even if thepressure is changing rapidly, and this reduces the number of stepsrequired in a simulation.

The invention may be applied to a “black oil” calculation or to acomposition calculation, but a streamline-based approach to solving a“black oil” calculation is described in a publication, SPE 51904,entitled “A Streamline Based Approach to Solution of Three-Phase Flow”,by L. Ingebrigtsen, et al, Feb. 14, 1999, published more than one yearbefore the priority date of this application. In a “black oil”calculation, only water, oil and gas are considered, and the individualhydrocarbon components are not considered, as they are in a full fluidcomposition calculation.

In one embodiment, a method of the invention is carried out using aStreamline Simulator (CS) and a finite difference simulator (FD). Thefinite difference simulator is used to provide the fully implicitsolutions along the streamlines. One suitable finite differencesimulator is the Schlumberger GeoQuest Eclipse 300 simulator, which is acommercial reservoir simulator. A suitable Streamline Simulator is theSchlumberger GeoQuest FrontSim simulator.

FIG. 3( a) is a flow diagram setting out the principal steps of a methodof the present invention. The method will be described with reference toa hydrocarbon reservoir but, in principle, the invention may beperformed on any volume that contains two or more fluid components.

Initially the CS simulator is initialised by loading user-defined datasuch as the pressure P, the water saturation S_(w), and the molefractions m_(t) of the various hydrocarbon components. Alternatively, ifthe gas saturation S_(g) and the mole fractions x_(i,j) of the j phasesof the i components are available, these can be loaded into the CSsimulator. The CS simulator contains a global grid of cells that areused in the calculations.

The FD simulator is then initialised for 1-D simulations using the fluidand component properties and relative permeability from the user data,and dummy data for the grid parameters and scheduling data. These dummyvalues are updated for each streamline during the simulation.

The CS simulator then instructs the FD simulator to calculate theinitial gas saturation, the initial pressure field, and the initialphase compositions for each compositional model. These calculations aredone in the flash calculation module of the FD simulator.

The FD simulator also calculates pressure derivatives.

The pressure field for a volume that defines the reservoir is thencalculated by the solver of the CS simulator. This is done using anIMPES-type solver, which has been modified to calculate the partialpressure field for each component of the reservoir fluid. Thepressure-dependent variables needed by the Jacobian are supplied by theflash module of the FD simulator

The streamlines are then computed from the pressure field.

The data that the FD simulator need to simulate one time step is thenprepared. The FD simulator has a fixed 1-D grid, and the streamlines arebundled up so that they fit into the fixed grid, as shown schematicallyin FIG. 4. The left hand of FIG. 4 is intended to represent streamlinethat extend from an injector to a producer as they would appear in theCS simulator (streamlines are generally three-dimensional, but are shownas two-dimensional in FIG. 4 for convenience). In order to solve alongthe streamline in the FD simulator, the streamline is mapped onto aone-dimensional grid, as shown on the right of FIG. 4.

The initial values for the FD simulator, such as the transmissibilityand depth of the grid, and the pressure, saturation and composition ofthe fluids, are set by the CS simulator. The pressure is computed bothglobally for each time step and locally along each streamline when thefluids are moved. [The positions of the producer(s) and injector(s) arealso defined.

The FD simulator then solves the pressure equation for one global timestep for its fixed 1-D grid which, as explained above, represents a setof streamlines. The FD simulator then passes the results of thissimulation to the CS simulator so that the values of the solutionvariables, the molar densities, the pressure and the saturations for theglobal grid cells can be updated. During this update, volumetricaverages of the values calculated on the streamlines are mapped back tothe global grid in the CS simulator.

FIG. 3( b) is a more detailed flow diagram showing one embodiment of theinvention.

First, the simulation is initialised. This initialisation includesinitialising two separate simulators, the CS simulator and the FDsimulator.

At step 10 user input data is read by the CS simulator. At step 11 theCS simulator prepares data to be input to the FD simulator, and at step12 this data is read into the FD simulator.

At step 13 the FD simulator is initialised for 1D simulations, typicallya small 1D grid. The size of the FD simulator grid does not depend onthe size of the CS simulator grid since the FD simulator is used tosolve for grids representing streamlines only. The fluid and componentproperties as well as the relative permeability tables are initialisedfor the FD simulator using data input by the user into the CS simulatorand passed to the FD simulator by the CS simulator. The rest of theinput data to the FD simulator (grid parameters, wells and productionschedule) is a set of dummy values for the initialisation. During thesimulation these variables will be dynamically updated for eachstreamline.

The CS simulator is then initialised at step 14. The initial gassaturation, pressure and phase compositions for the CS simulator arecomputed using the FD simulator flash calculation module at step 15, andare supplied to the CS simulator. The FD simulator flash module alsosupplies the CS simulator with data for the phase mobility, density, andcompressibility as well as the derivatives of these quantities withrespect to pressure.

Next the CS simulator computes the pressure field within the reservoirat step 16. The pressure solver in the CS simulator is an IMPES typesolver and few modifications have been made for the compositional case.The pressure dependent variables needed for the matrix are supplied bythe FD flash module. From the pressure field the fluxes are determinedand streamlines are computed at step 17. The streamlines may be computedin any suitable manner, for example using the method described above.

In order to use the FD simulator to solve for the compositions andpressure, the streamlines are broken down into 1D grids defined bytransmissibility, pore volume and depths for the grid cells. This isdone in the CS simulator, as step 18. The gravity lines are also brokendown into 1D grids, again at step 18.

Each of these 1D grids can be solved individually by the FD simulator.For performance reasons, the CS simulator will fill the FD grid with aset of streamline grids. An injector is positioned in the first cell anda producer in the last cell of each part of the grid representing asingle streamline, as shown schematically in FIG. 4.

At step 19, the FD simulator solves one time step for the set ofstreamlines in the FD grid. The results of this step are sent to the CSsimulator. The resulting compositions and pressure determined at step 19are mapped onto the CS grid at step 20. This procedure is repeated untilall streamlines have been solved.

Once all streamlines have been solved and mapped onto the CS simulatorgrid, the FD simulator performs, at step 21, a flash calculation on theresulting averages in the CS simulator grid. As for the initial flashcalculation at step 15, the FD simulator flash will supply the CSsimulator with values for the phase mobility, density andcompressibility as well as the derivatives of these quantities. However,the values calculated at step 21 are updated values, at the end of atime step. These values may be output for review by an operator, if thisis desired.

At step 22 there is a determination whether the simulation is complete.If the result of this determination is “no”, a time counter isincremented at step 23, and the pressure is then re-calculated at step16. The pressure is recalculated from the updated values computed by theFD simulator at step 21.

Steps 16 to 23 are then repeated until the simulation is complete and a“yes” determination is obtained at step 22, whereupon the process endsat step 24. The detailed mathematical basis for the present inventionwill now be described.

The gravity lines are handled very much in the same way as thestreamlines. They are mapped to the FD simulator, and results areobtained in the same way as for the streamlines. The main difference isthat there are no wells defined in the gravity lines since only thegravity segregation of the fluids need be taken into account. Theinitial condition for the gravity lines are the results of thestreamline simulator.

From equation (1), the component conservation equations for eachcomponent i in the FD simulator can be written in terms of n_(c)+2primary solution variables for pressure and component molar densitiesX=(p, mi, mw), where i=1 . . . n_(c), and mw is the molar density ofwater.

$\begin{matrix}{{{\frac{\partial}{\partial t}m_{i}} + {\nabla{\cdot ( {\sum\limits_{{j*g},w,o}{x_{i,j}{\lambda_{j}\lbrack {{\nabla p_{j}} - {\rho_{j}\overset{\rho}{g}}} \rbrack}}} )}} + q_{i}} = 0} & (12)\end{matrix}$for i=1 . . . n_(c), where j denotes the phase and λ_(j) is the mobilityof the j^(th) phase,

$\lambda_{j} = \frac{k_{r\; j}b_{j}}{\mu_{j}}$(with k_(rj) the relative permeability of phase j). A similar equationcan be written for water.

The finite volume discretisation for component and water molarconservation equations then lead to n_(c)+1 equations of the form

$\begin{matrix}{{\frac{( {( {V_{p}m_{t}} )^{n + 1} - ( {V_{p}m_{t}} )^{n}} )}{\Delta\; t} + {\sum\limits_{j}F_{i,j}} + q_{i}} = 0} & (13)\end{matrix}$

In streamline mode, the ID spatial discretisation in equation (13) is intime of flight (TOF) along the streamlines. If the total number ofhydrocarbon moles per unit pore volume is

$m_{t} = {\sum\limits_{i = 1}^{n_{c}}\; m_{1}}$then the phase saturations for oil, water and gas, S_(w), S_(w), S_(g),are defined in terms of the liquid and vapour mole fractions L, V andthe phase molar densities of oil, water and gas, b_(w), b_(o), b_(g), as

$\begin{matrix}{{S_{w} = {m_{w}/b_{w}}},{S_{o} = {{Lm}_{t}/b_{o}}},{S_{g} = {{Vm}_{t}/b_{g}}}} & (14)\end{matrix}$where m_(t) is the total hydrocarbon molar density, and so the phasevolume per unit pore volume is given by

$\begin{matrix}{u_{T} = {\frac{m_{w}}{b_{w}} + \frac{{Lm}_{t}}{b_{o}} + \frac{{Vm}_{t}}{b_{g}}}} & (15)\end{matrix}$

The FD simulator solves the additional volume balance constraintV_(p)=V_(f)=V_(p) u_(T), where V_(p) is the pore volume, V_(f) is thefluid volume, and u_(T) is the phase volume per unit pore volume, andtotal compressibility of the fluid is given by

$\begin{matrix}{C_{t} = {{- \frac{1}{u_{T}}}\frac{\mathbb{d}u_{T}}{\mathbb{d}P}}} & (16)\end{matrix}$

This is used to construct the pressure equation (3). The FD simulatorflash package also supplies values and derivatives for other extensiveproperties in equation (3) such as μ_(j), the viscosity of phase j, andρ_(j).

The CS simulation is initialised by calling the compositionalequilibrator, and is flashed onto the global grid. The initial dataconsists of the pressure P, the saturation of water S_(w) and molefractions z_(i), or both S_(w) and the gas saturation S_(g) and molefractions x_(i,j) if available. Since the simulator works with specificmolar density primary variables m_(i) it is important to predict m_(i)in order to preserve material balance when the solution is interpolatedback and forth between the global grids and streamlines. After solving asimulator time step, the solutions are then pulled back onto the globalgrid and reflashed before computing the new streamlines. Materialbalance errors and the and the work involved in this reflash isminimised by initialising the FD streamline solver and global flash withsaturations S_(o), S_(g) and S_(w) and liquid and vapour compositionsx_(i,o) and x_(i,g), so the vapour fraction V=S_(g) b_(g)/m_(i) ispredicted. The initial molar phase densities are computed asb_(j)=b_(j)(x_(i,j),P) and b_(w)=b_(w)(P) so that the component molardensities are given bym _(i) =x _(i,o) S _(o) b _(o) +x _(i,g) S _(g) b _(g)  (17)

The interpolation of saturations S_(j) and mole fractions x_(i,j)between the streamline and global grids attempts to conserve the molarmasses x_(i,j) S_(j) b_(j) of components in each phase as well as themolar mass of water m_(w) while pressures are pore volume weightedaverages. The material balance errors are then minimal and uT should beapproximately 1. If uT is significantly different from 1, the predictedS_(o), S_(w), S_(g) in equation (15) obtained from the FD simulatorflash package are normalised by uT, leading to material balance error inthe streamline simulation.

The compositional simulator runs efficiently since the streamlines areone-dimensional. Streamlines are concatenated into a I-D list andpresented to the simulator as a bundle so that key parts of the FD canreadily be optimised for a 1-D solution. Each streamline in the bundlemay be converged separately; the convergence criteria in the simulatorare based on achieving tolerances on solution changes in p and mi.

In the examples described in detail above, the streamlines have beensolved fully implicitly. Alternatively, the streamlines may be solvedusing adaptive implicit method (AIM) solutions along streamlines.

Examples illustrating the improvement provided by the invention will nowbe described.

1. SPE9 Simple Model

The 9th SPE Comparison Problem, as presented by J. Killough at the 13thSymposium of Reservoir Simulation, San Antonio, Tex. (1995) wasconverted from black oil to a compositional problem. The problem relatesto heterogeneous reservoir with a single water injector and 24 oilproducers whose operation was converted to reservoir volume control. Inthe simulation a 9000 cell grid (24×25×15) is used. The producers arecompleted (open) in the layers 2,3,4, and the injector was completed inlayers 11 to 15. The reservoir is initialised with an oil water contactin the reservoir. The simulation is run for a period of 900 days. ThePVT data for 6 components, having from one to 20 carbon atoms (C1–C20),is used with single stage flash to standard conditions of 60° F. and14.7 psi.

FIG. 5 shows the field gas production rate (FGPR) as simulated using aCS model of the present invention and an FD model for comparison. Itwill be seen that the material balance error predicted by the CS methodof the invention is less than 0.6%.

2. SPE Refined Model

The SPE9 comparison model was refined to a 72×75×45 grid with the same 6component PVT data as for the SPE simple model described above to give a243,000 cell model. For models of this size, the streamline approach maybe much faster than a finite difference simulation. A simplifiedproduction schedule was employed with the production wells under aconstant reservoir volume control of 1500 Stb/day and the water injectormaintaining a bottom hole pressure (BHP) of 4500 psi over a period of1000 days. The producers were completed (open) in layers 4 to 12). FIG.6 compares field oil production rates (FOPR), field gas production rates(FGPR) and field water production rates (FWPR) for a CS simulation andan FD simulation. It can be seen that the oil recovery predicted by theCS simulation is slightly lower than the oil recovery predicted by theFD simulation.

FIG. 13 shows the memory capacity and CPU time required to perform asimulation for this model, using a CS simulation of the invention andusing a conventional FD simulation. The CS simulation used a 243,000cell model, and extended over a simulation period of 1,000 days. ThisFigures shows that a CS simulation of the invention is quicker and needsless memory that a conventional FD simulation. The CS simulationrequired 14 hours of CPU time (on an IBM SP2) and 330 MB RAM, whereasthe conventional FD simulation required 75 hours CPU time and 650 MBRAM. Thus, even with no tuning, a CS simulation of the present inventionruns in less than a fifth of the CPU time required for an FD simulation,and needs approximately half the memory capacity.

FIGS. 7 and 8 illustrate the saturation (FIG. 7) and the streamlines(FIG. 8) for the top layer of the reservoir simulation for the SPE 9refined model after 500 days. In FIG. 7 the saturation is shaded todenote oil, gas or water. As indicated in the Figure, a region near theinjector extending across the entire breadth of the top layer of thereservoir (that is, into the paper), and over approximately 0.2 of thewidth of the top layer of the reservoir (that is, across the paper)contains predominantly water, The remainder of the top layer of thereservoir contains predominantly oil, except for a narrow region at theboundary of the water-containing region.

The streamlines in FIG. 8 are shaded to indicate gas saturation. whereasthe is where the streamlines are shaded to indicate gas saturation. Thegas saturation is close to zero near the injector, and increases awayfrom the injector, generally to around 0.1 although it reaches o.2 or0.3 in some regions. The material balance error is 0.03% for the CSsimulation.

3. The PUNQ Simple Model

The PUNQ-Simple model of F. J. T Floris et al, European Conference onthe Mathematics of Oil Recovery (1998) was converted to a 5-componentcompositional model and initialised with an oil water contact in thereservoir. The coarse simulation model has 1761 active cells and the rimaquifer was removed for the purpose of comparing the finite volumesimulation and the streamline simulation. The reservoir is depletingwith the wells put under reservoir volume rate controls for a 5000 dayperiod from January 1967. The well separator conditions assume a singlestage flash to standard conditions. FIG. 9 shows the field oilproduction rates (FOPR) predicted by a conventional FD simulation and bya CS simulation of the invention. It can be seen that CS simulationpredicts an early decline in production when the BHP constraint isachieved. (Initially in the simulation the produces are controlled tohave a constant rate, and a limit for the maximum bottom hole (BHP) isset. The BHP increases during the simulation, and will eventually reachthe set limit. At this point the target for the well is changed tomaintain a constant BHP, and the production rate decreases.

FIGS. 10 and 11 show the saturation plot (FIG. 10) and the streamlineplot (FIG. 11) for the CS model after 5000 days when free gas hasappeared in the reservoir. The indicated region in the saturation plotof FIG. 10 contains predominantly oil, and elsewhere there as anoil/water mix (with the water content generally increasing away from thepredominantly-oil region. The gas saturation in FIG. 11 wasapproximately 0.2 in a region that corresponds generally to thepredominantly-oil region of FIG. 10, although it reached 0.25–0.3 nearproducers Pro 1, Pro 4 and Pro 5. Elsewhere, the gas saturation waslower, generally in the range 0.05 to 0.15.

4. PUNQ Complex Model

As an illustration of an application of the invention to a geologicalscale model that it is impractical to solve using a conventionalsimulation without first upscaling, the million cell PUNQ complex model(the SPE 10^(th) comparative solution project on upscaling) was run. Inthis model the wells are arranged in a 5 spot pattern. A centralinjector injects water at a reservoir volume rate of 4000 Stb/day andfour producers operate under bottom hole pressure constraints of 4000psi. The reservoir is located between 12,000 and 12,170 feet and isinitialised with a datum pressure of 6000 psi at 1200 feet with an OWCat 2000 feet. The six component fluid described in Example 1 is used.FIG. 12 illustrates the streamlines after 5110 days. The streamlines areshaded to indicate the C10 molar density, and it can be seen that thisincreases from approximately 0.25 at the centre of the layer toapproximately 1 at each of the producers P1 to P4

The invention has been described above with reference to separate CS andFD simulators. In practice however, the invention may conveniently becarried out by embodying the CS simulator and the FD simulator in asingle program that is run on any suitable processor.

FIG. 14 is a schematic block diagram of an apparatus 25 according to thepresent invention. The apparatus comprises an input device 26 forreceiving input data to be processed and for supplying to the remainderof the apparatus data t the appropriate time for processing.

The apparatus further comprises a streamline simulator 27 and an FDsimulator 28. The FD simulator may be a fully-implicit FD simulator, orit may be an adaptive implicit FD simulator.

Finally the apparatus comprises an output device 29, such as a printer,a visual display unit, or a memory.

FIG. 15 illustrates a programmable system embodying the apparatus ofFIG. 14, for performing a method according to the present invention,such as the method illustrated in FIG. 3( a) or FIG. 3( b). The systemcomprises a programmable data processor 30 with a program memory 31, forinstance in the form of a read only memory ROM, storing a program forcontrolling the data processor 30 to perform, for example, the methodillustrated in FIG. 3( a) or FIG. 3( b). The system further comprisesnon-volatile read/write memory 32 for storing, for example, any datawhich must be retained in the absence of power supply. A “working” or“scratchpad” memory for the data processor is provided by a randomaccess memory (RAM) 33. An input interface 34 is provided, for instancefor receiving commands and data. An output interface 35 is provided, forinstance for displaying information relating to the progress and resultof the method. User-define data may be supplied via the input interface34 or may optionally be provided by a machine-readable store 36.

The program for operating the system and for performing the methoddescribed hereinbefore is stored in the program memory 31, which may beembodied as a semi-conductor memory, for instance of the well-known ROMtype. However, the program may be stored in any other suitable storagemedium, such as magnetic data carrier 31 a (such as a “floppy disc”) orCD-ROM 31 b.

1. A method of determining fluid flow in a volume containing three ormore hydrocarbon components, the method comprising: determining apressure field for the volume; gridding the volume; determining one ormore streamlines from the pressure field; bundling the, or each,streamline into a one dimensional grid; and solving for, and storing, afull fluid composition along the, or each, streamline.
 2. A method asclaimed in claim 1 and further comprising solving for the pressure alongthe, or each, streamline.
 3. A method as claimed in claim 2 furthercomprising the steps of: determining an updated pressure field for thevolume; determining one or more streamlines from the updated pressurefield; and solving for a full fluid composition of the reservoir alongthe, or each, streamline determined from the updated pressure field. 4.A method as claimed in claim 3 wherein the step of solving along the, oreach, streamline comprises using a finite difference technique.
 5. Amethod as claimed in claim 3 wherein the step of solving along the, oreach, streamline comprises using a fully implicit finite differencetechnique.
 6. A method as claimed in claim 3 wherein the step of solvingalong the, or each, streamline comprises using an adaptive finitedifference technique.
 7. A method as claimed in claim 2 furthercomprising the steps of: determining an updated pressure field for thevolume; determining one or more streamlines from the updated pressurefield; and solving for a full fluid composition and pressure along the,or each, streamline determined from the updated pressure field.
 8. Amethod of processing data relating to a hydrocarbon reservoir containingthree or more hydrocarbon components, the method comprising the stepsof: determining a pressure field for a volume indicative of thereservoir; gridding the volume; determining one or more streamlines fromthe pressure field; bundling the, or each, streamline into a onedimensional grid; and solving for, and storing, a full fluid compositionof the reservoir along the, or each, streamline.
 9. A method as claimedin 8 further comprising solving for the pressure of the reservoir alongthe, or each, streamline.
 10. A method as claimed in claim 9 furthercomprising the steps of: determining an updated pressure field for thevolume indicative of the reservoir; determining one or more streamlinesfrom the updated pressure field; and solving for the fluid compositionof the reservoir or for the fluid composition and the pressure along theor each streamline determined from the updated pressure field.
 11. Amethod of processing data relating to a hydrocarbon reservoir containingtwo or more hydrocarbon components, the method comprising the steps of:determining a pressure field for a volume indicative of the reservoir;determining one or more streamlines from the pressure field; and solvingfor the fluid composition of the reservoir using compositionalsimulation, along the, or each, streamline; solving for the pressure ofthe reservoir along the, or each, streamline; determining an updatedpressure field for the volume indicative of the reservoir; determiningone or more streamlines from the updated pressure field; solving for thefluid composition of the reservoir or for the fluid composition and thepressure along the, or each, streamline determined from the updatedpressure field; and further comprising the further step of controllingproduction of the reservoir on the basis of the results of the step ofsolving along the, or each, streamline.
 12. A method as claimed in claim11 wherein the step of solving along the, or each, streamline comprisesusing a finite difference technique.
 13. A method as claimed in claim 11wherein the step of solving along the, or each, streamline comprisesusing a fully implicit finite difference technique.
 14. A method asclaimed in claim 11 wherein the step of solving along the, or each,streamline comprises using an adaptive finite difference technique. 15.An apparatus for determining fluid flow in a volume containing three ormore hydrocarbon components, the apparatus comprising: means fordetermining a pressure field for the volume; means for gridding thevolume; means for determining one or more streamlines from the pressurefield; means for bundling the, or each, streamline into a onedimensional grid; and means for solving for a full fluid compositionalong the or each streamline.
 16. An apparatus for processing datarelating to a hydrocarbon reservoir containing three or more hydrocarboncomponents, the apparatus comprising: means for determining a pressurefield for a volume indicative of the reservoir; means for gridding thevolume; means for determining one or more streamlines from the pressurefield; means for bundling the, or each, streamline into a onedimensional grid; and means for solving for a full fluid composition ofthe reservoir along the, or each, streamline.
 17. An apparatus asclaimed in claim 16 wherein the means for solving the fluid compositionalong the, or each, streamline is a finite difference simulator.
 18. Anapparatus as claimed in claim 17 wherein the finite difference simulatoris further adapted to solve for the pressure along the, or each,streamline.
 19. An apparatus as claimed in claim 18 wherein the finitedifference simulator is a fully implicit finite difference simulator.20. An apparatus as claimed in claim 18 wherein the finite differencesimulator is an adaptive implicit finite difference simulator.
 21. Anapparatus as claimed in claim 20 wherein the means for determining thepressure field is a streamline solver.
 22. An apparatus as claimed inclaim 21 comprising a programmable data processor.